Random matrix ensembles: Wang-Landau algorithm for spectral densities 
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We propose a method based on the Wang-Landau algorithm to numericaUy generate the spectral 
densities of random matrix ensembles. The method employs Dyson's log-gas formalism for random 
matrix eigenvalues and also enables one to simultaneously investigate the thermodynamic properties. 
This approach is a powerful alternative to the conventionally used Monte Carlo simulations based 
, on the Boltzmann sampling, and is ideally suited for investigating /3-ensembles. 
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I. INTRODUCTION 



Random matrix theory (RMT) was introduced in physics by Wigner to describe the spectral properties of heavy 
nuclei. It is now a huge subject with applications in areas as varied as mesoscopic physics, quantum field theory, 
number theory, wireless communications and econophysics [iHEl- 
D ' The statistics of eigenvalues plays a crucial role in the study of random matrix ensembles (RME). It not only 
^ I encodes valuable information about a given ensemble but also enables one to calculate important observables related 
^ I ^ ■ to the eigenvalues. Quite often one needs to evaluate the eigenvalue densities numerically. This is because the 
' analytical expressions for the densities are available only for certain special ensembles and/or in asymptotic limits. In 
^ . addition, the numerical simulation of spectra facilitates the verification of known analytical results. The conventional 
^ ; ' numerical schemes for generating eignvalue densities rely on either the diagonalisation approach using the matrix 
J model whenever available, or implementing the Boltzmann-sampling-based Monte Carlo simulation 6]. The latter 
^ ' uses Dyson's log-gas picture for the eigenvalues of RME [l|, [?[ . We propose here a method based on Wang-Landau 
I ] algorithm (WLA) Q which serves as a powerful alternative to these approaches. 
' ^ ■ Our proposed method relies on calculating the microcanonical density of states of the Dyson log gas associated with 
the given random matrix ensemble. The density of eigenvalues or distribution of any related observable (function of 
eigenvalues) can be obtained from this information. RME are characterised by the Dyson index /3 which plays the role 
of inverse temperature in the log-gas thermodynamics. The superiority of WLA over the conventional schemes lies 
in that one can evaluate the spectral density as well as the distributions of related observables for multiple /? values 
in a single simulation. The method, therefore, can be viewed as a natural tool in exploring the /3-ensembles [9l-[ll|. 
Moreover, WLA provides the additional advantage of producing density of states of the log gas during the simulation, 
which can be used to investigate its thermodynamic properties [l|, 0| . 



in 

II. WANG-LANDAU ALGORITHM 

o 

[ Wang-Landau algorithm [8| is a non-Boltzmann sampling method which was originally applied to standard statistical 
systems like Ising and Potts models which exhibit discrete energy states. WLA aims at estimating accurately the 
microcanonical density of states g{E) (up to a multiplicative constant) by performing a random walk in the energy 
k>l space. The key idea behind this algorithm is that if during the random walk the configurations x with associated energy 
r> ■ E are sampled with a probability proportional to l/g{E), the corresponding histogram H{E) of visited configurations 
\ becomes flat. The density of states being unknown at the beginning of the simulation, one starts with an initial guess 
" " " for g{E) which is converged towards the true density in an iterative manner. Since the method gives direct access to 
the density of states of the system, which is independent of the temperature, one can calculate various thermodynamic 
averages by canonical reweighting at arbitrary nonzero temperature Q. There has been several improvements and 
refinements to the WLA by various authors [l3 - [l7t , consequently it has been successfully applied to systems possessing 



continuous energy spectra also, e.g., liquid crystals, polymers, biomolecules etc. See for example |15| . where the 



authors propose a hybrid algorithm based on the Wang-Landau and transition matrix Monte Carlo methods. In our 
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FIG. 1: (a) Log of the unnormalised density of states (b) Density of eigenphases for Dyson ensemble, eq. ([Sjl with w(e'^) = 1 
(c) Density of observable A = Sj/N. Uniform distribution pe{6) = (2tv)~^ is obtained in (b) for each /3 value. In (c) the 
symbols are the simulation results while the solid lines represent the analytical results given by eq. (|18[) . 



calculations we implement the variant [l6l.[T7j of WLA which facilitates, in general, a faster convergence rate than 
the conventional scheme. 

From the information of g{E) we can immediately evaluate the canonical probability at arbitrary nonzero temper- 
ature T = l/(fcB/3) as 

ks being the Boltzmann constant. The sum in the above equation runs over the sampled energy points in the energy 
window of interest. The canonical average of an observable depending explicitly on energy, say Q{E), can then be 
easily calculated as 

{Q) =Y.Q{E)Pp,E{E). (2) 

E 

For example, the moments of energy (i?") can be evaluated as J2e where n = 0,1,2,.... By definition 

we have (S") = 1. The average energy is given by the first moment, U — (E^). The specific heat, which is related to 
the second cumulant of energy, is obtained using Cy = kB(3^ ((^^) ~ i^^)^)- 

For evaluation of the quantities which do not depend explicitly on E, say A, there are two ways one can proceed. 
The first approach is based on an extension, where instead of performing the random walk in only .E-space, one 
performs the random walk in two-dimensional {E, A) space with probability proportional to 1/Q{E, A), G{E, A) being 
the two-dimensional density of states . The distribution of observable A can be calculated at the end of simulation 
using 

which can be further used to calculate the average (A) using 'Ej^APfj_E{A). However, performing the random walk 
in two-dimensional space is usually computationally expensive, being plagued with convergence issues etc. [Tsj . The 
alternative approach relies on estimating g{E) first and then performing a production run [Tol. [20j. The production 
run comprises a post-simulation employing "Wang-Landau resampling" using the already converged g{E). During 
this run, a histogram H{E,A) of the visits in the {E,A) space is generated. The joint density Q{E,A) can then be 
obtained using [1^] 

g{E,A)^g{E)n{E,A). (4) 



Eq. Q can finally be used to study the desired observable. 
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FIG. 2: Density of eigenphases for wie'") = |1 + 2e'® + 3e'^®|"^'' and N = 10. Solid lines are the result of Wand-Landau 
simulation. For comparison we also show the Boltzmann-sampling based Monte Carlo results (symbols). WLA gives the 
densities for arbitrary /3 > in a single simulation whereas the conventional Monte Carlo needs separate simulation for each /3. 

III. DYSON LOG GAS 

Dyson in his pioneering work provided a physical picture to the eigenvalues of RME [7^. He showed that the 
eigenvalues can be identified with the positions of charged particles interacting via two-dimensional Coulombic repul- 
sion and executing Brownian motion. In addition a position dependent potential may be present. Dyson's log-ga s 
model P, 0, 01 is of fundamental importance to the field of RMT and has been applied in various contexts [2l'-'25'] . 
Under this formalism the three invariant ensembles of random matrices, viz. orthogonal, unitary and symplectic [1], 
turn out to be the equilibrium states of this system at three special temperatures (3~^ = 1, 1/2, 1/4, ks being set equal 
to 1. Dumitriu and Edelman, by dropping the invariance requirement, proposed matrix models which correspond to 
general /3 > 0. Moreover, very recently it has been shown that by considering a sp ecial diffusive matrix model, 
invariant /3-ensembles (/3 restricted to certain continuous domains) may be realised [lOl [Tl|. 

In RMT one usually deals with the ensembles of unitary matrices and Hermitian matrices, which share the following 
generic structure for the joint probability density (JPD) of eigenvalues {xjlJ — 1, .■.,-/V): 

Pxilx}) = n IX. -XkfU w{xi)- (5) 
]<k I 

Here {%} = {xi: Xn} and Z = Y[j<k \Xj ~ Xk\^ Yii ^{xi)dxi is the partition function, R being the appropriate 
domain of integration. \xj — Xk\^ is the characteristic random matrix eigenvalue-repulsion term, and 'w{x) is the 
weight function which decides specific features of a given ensemble. Note that, for brevity, we use the term Hermitian 
to refer to real symmetric and self-dual quaternion matrices also. 

In the unitary case the eigenvalues are located on a unit circle in the complex plane, i.e., they are of the form 
e'^^, dj € [— 7r,7r] being the eigenangles (or eigenphases). We will use x below to refer interchangeably to e** and 6. 
w(e*^) = 1 corresponds to the important case of Dyson's (uniform) circular ensemble In the Hermitian case Xj = 
are positions on the real line. Jacobi family of RME is obtained when w{x) is chosen as weight functions associated 
with the classical orthogonal polynomials. For example, w{x) = e~^^ 2./3(a-i-i)/2g-/3a;/2^ (X_2;)^(a-i-i)/2^X-|-a;)'^(''+^)/^ 
lead to the Gaussian, Laguerre and Jacobi ensembles of random matrices [l|, |3, d, [13] . In these three cases the xj lie 
in R = (—00,00), [0,00) and [—1, 1] respectively. The restrictions on a,b are determined from the normalisability of 
the weight functions. 

The JPD of eigenvalues can be used to extract important information about the ensemble. For example, we can 
calculate the (marginal) density of eigenvalues, viz., Pxix) = N~^{J2j Hx — Xj))K- Here ( )h represents the ensemble 
average with respect to the JPD in eq. ^ . (x) provides the information about average behaviour of the eigenvalues 
and therefore reveals the general features of a given ensemble. Similarly, the distribution of an observable A — A{{x}) 
depending on x's can be obtained using /a(^) = {5{A — A{{x})}r- 
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The JPD of eigenvalues in eq. ([5]) can be recast as a canonical Boltzmann-Gibbs probability density: 

P^{{x})^Z-'e-f'''^^^^\ (6) 
Here E represents potential energy of the system, 

i?=-^ln|x,-Xfe|+E^fe)- (7) 

j<k j 

In the log-gas picture the first term in the above equation represents the two-dimensional Coulombic repulsion between 
the charges, while V{x) represents the local binding potential seen by the charges [H, i, 0, [H HI . V{x) is related to 
the weight function w{x) in eq. ([S]) as 

V{x) = -r'^riw{x). (8) 

Eq. ([7]) provides the energy function which we use in WLA to perform the random walk in i?-space and to derive the 
density of states g{E). g{E) can eventually be used to generate the density of eigenvalues and distributions of related 
observables. We give the details in the next section. 

It can be shown that for the above mentioned classical weight functions, E is bounded from below, i.e., there 
exists an Eq such that E > Eq for all Xj El- For the Dyson's circular ensemble, minimum energy configuration 
corresponds to the case when the angular separation between any two adjacent charges on the unit circle is equal. In 
this case we have [l[ 

£;o = -^A^lniV. (9) 

For the Hermitian case (eigenvalues as positions on real line) the minimum energy configuration {x^i \ x"^} coin- 
cides with the zeros of classical orthogonal polynomials. For the weight functions given above, x^^^'s are precisely the 
zeros of Hermite: Hn^sc), associated Laguerre: L^^\x), and Jacobi: Pj^'^\x) polynomials P, HI]. The corresponding 
expressions for Eq are [l|, 



,10, 

for the Gaussian case, 



2 



2 

3 



,. , a+1, T{N + a+l) 

-> lii(]+a) In— (11) 

^ 2 ^ 2 r(a -M ^ ' 

3 

for the Laguerre case, and 

i.o^-^(^ + - + ^+^)ln2-(A^-l)ln^(^+^)^^^ + " + ^+^) 
2 ^ ' T{2N + a + h+l) 

3 3 3 

a + 1 T{N + a + l)r(iV + a + b + I) ^ + I T{N + h+ l)T{N + a + b+l) 



r(a+l)r(27V + a + 6+l) 2 r(5+l)r(2Ar- 



(12) 



for the Jacobi case. 

In the above examples for the classical ensembles we have chosen w(x) in a way that V{x) in eq. (|5]) is /3-independent. 
However, often this is not the case [1, H, [13, HH] • For instance, the density of transmission eigenvalues, which govern 
the statistics of conductance (or transmission) in chaotic cavities, is given by [1, 

Pt{{T}) oc n \T,-n\f^\{T^^''-''+'^'^-\ (13) 
j<k I 
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FIG. 3: Density of eigenvalues for (a) Gaussian: = 10, (b) Laguerre: A'^ = 5, a = 2, and (c) Jacobi: N = 9,a — 1,6 = 3 
ensembles (Symbols: Simulation; Solid lines: Analytical [2J| ). The legend applies to all three boxes. 

Here < T,- < 1, and AI,N{M > N) represent the number of modes in the two ideal leads connected to the 
cavity. V{T) in this case is —((Af — + l)/2 — InT, which makes the energy £^ in eq. ([7]) /3-dependent too. Asa 
consequence, for each /3 we will need to do a separate simulation. This completely nullifies the usefulness of WLA over 
the conventional Monte Carlo scheme. However, we can circumvent this problem by considering the transformation 
Tj = e~*\0 < tj < oo. The factor coming from the Jacobian of transformation gets rid of the unwanted —1 in the 
exponent and we obtain 



]<k 



-tk\l3 



-f}(M-N+l)ti/2 



(14) 



WLA can now be implemented using the /3-independent energy function, eq. ([7]) with Xj = e *5 and V — {b+l)tj /2, b 



M ~ N . Moreover, in this case 

Eq = -{N - l)ln 



T{N + l)T{N + h) 



T(2N + h) 



E 



AT _ j .. . 

Y^—^\n{N + j+b-l)-Y,'-—Hj-l 



j-2N + 2 



Inj 



j 



iln(j- + 6)-- 



1 ^^ TjN + b+l)r{N + b) 
T{b+l)T{2N + b) 



(15) 



The term — ln(j — 1) in the sum above should be taken (in a limiting sense) as for j = I. With the /3-independent 
energy function we can obtain the density pt{t) for multiple values of /3 in a single simulation. The density of T's can 
then be trivially obtained using the relation pt{T) = e^pt{t). 

As our next example we consider the JPD for /3-Wishart ensemble as derived in [T]| : 



^A({A})(xniA,-Afei^n^'^^^^ 



i . |(M-W+l-5)-(l-|) 



(16) 



Here < Aj < 00,(5 > and N/M < 1. It has been shown that for large M,N, a crossover between Marcenko- 
Pastur and the Gamma distribution (/? — 0) is observed for /3 ^ 1/M [llj. Similar to the previous example, 
the above JPD as such is not suitable for the implementation of WLA. We therefore carry out the transformation 

2 / 6 

Xj — Py-' ,0 < Uj < oo. This gives us the following equation for JPD of j/'s: 



Py{{y}) « n 



Vk 



i(M-N+l-5) 



(17) 



j<k 



This is essentially the Laguerre JPD in variables y"^^^ with a in the effective weight function 
(y2/<5^^(a+i)/2 exp(— /3y^/'^/2) identified as M — N — S. The corresponding i?-function is given by eq. ([7|) with x = vV^ 
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FIG. 4: (a) Density of transmission eigenvalues for M = 15, A*' = 4, (b) Conductance distribution for M = N = 2, and (c) 
Average conductance as a function of /3 for several M,N values. Symbols: Simulation; Solid lines: Analytical p5l |. 

and V = 2^^y^^^ — S^^{a + 1) Inj/j. The ground energy Eq follows from eq. dTTI) by setting a = M — N — S. The 

expression of density in original variables can be obtained as p\{X) — {2f3)~^5y^~^py{y). We would like to remark 
that one can also consider the case 5 = in eq. (fT6)) . for which the JPD corresponds to the classical Wishart ensemble. 
In this case we carry out the transformation Xj — j3e~y\ —oo < yj < oo, similar to the one considered in transmission 
eigenvalue JPD above. 

IV. IMPLEMENTATION OF THE ALGORITHM 

While investigation of the thermodynamic behaviour of these ensembles is interesting in itself [l], 0, HBl , our focus 
here is on the evaluation of density of eigenvalues. 

A successful implementation of WLA using the energy function, eq. ([7]) relies primarily on a reasonable choice of 
the {E,x) space to be explored, i.e., the range (windows) for eigenvalues/eigenangles and energy. For the energy 
window (We), it is advantageous to shift in eq. ^ as E ^ E — Eq if Eq is a priori known. Thus we perform the 
simulation with E — Eq as the energy function with the lower end of We set to zero. The energy being unbounded 
from above, fixing a reasonable upper cut off for the We requires some hit and trial. We can get some idea from a 
pre Monte Carlo run where energies are generated for random configurations. This is also helpful when there is no 
a-priori information about Eq, and one has to set a reasonable lower cut off also. As an example of such a case we 
have considered the weight function w{e^^) = |1 + 2e'^ + 3e'^^|~^'^ in eq. ([5|). Another option is to implement the 
self-adaptive method suggested in (isj . 

The X" window (W^) can be trivially fixed for ensembles which have a finite support for the eigenvalues/eigenangles, 
e.g., Jacobi {xj e [—1,1]) and circular {9j G [— tt, tt]) ensembles. However, when the eigenvalues do not have a finite 
support, e.g. in Gaussian and Laguerre ensembles, we have to carefully choose the window. The choice of W-,^ in 
these cases should be such that the g{E) in the We does not change appreciably with further increase in the W^ 
width. One can also get a crude estimate using the large TV-asymptotic result, if available. For example, in the case 
of Gaussian weight e~'^^ in eq. ([S]), the eigenvalues give rise to Wigner semicircle of radius V2N in the large N 
limit. Therefore we can restrict the W^ size to be of this order. Afterwards, the windows are divided into reasonable 
number of bins for sampling in the (E, x) space. In general, too many bins will hinder the convergence, while too 
little bins will lead to poor results. 

As mentioned before, we implement the variant of WLA for the simulation [iB, Once the g{E) has 

converged to the desired accuracy Q, we perform the production run to obtain the joint histogram 'H{E,x)- Owing 
to the symmetry of x's we update 'H{E,x) by simultaneously tracking all the eigenvalues instead of just one. This 
speeds up building of the histogram. The density of x's can finally be obtained using eqs. ([3]) and ([4]). We can also 
obtain the distributions or averages of multiple observables which depend on the eigenvalues by maintaining joint 
histogram for each of them; as has been done in some of the examples considered. 
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FIG. 5: Density of eigenvalues for /3-Wishart ensemble, eq. ([16]), for M — 100, A'^ = 50 and S = 1 (Symbols: Simulation, Solid 
lines: Analytical [lll|). A crossover from Gamma distribution to Marcenko-Pastur density is observed for j3 ~ 1/AI. 



RESULTS AND DISCUSSION 



The results for Dyson's circular ensemble, w(e*-^) = 1, are shown for TV = 2 in fig [U Figures l(a)-(c) show the 
logarithm of density of states g{E), density of eigenphases 9j, and distribution of observable A = 9j. The 

eigenangle density for each /3 value is a uniform distribution [l|. In fig. 1(c) along with the simulation results (symbols) 
we also show the N — 2 analytical results (solid lines): 



.fA{A) 



r 1 nA 

TT Z 

---l^\2A-sm2A\, 

1 A 
— (2 — cos A) cos^ — , 

TT 2 



A 



2 sin 2 A sin AA 



12 



= 
/3 = 2, 

/3 = 3, 

/3 = 4. 



(18) 



In fig. [21 the density of eigenphases is shown for the weig ht function w{e''^) = |1 + 2e** + Se*^^!-^'^ for N = 10. The 
solid lines represent the results from the Wang-Landau simulation. For comparison we also show the results of Monte 
Carlo simulation based on Boltzmann sampling 

Fig. [3| shows the eigenvalue densities for (a) Gaussian, (b) Laguerre and (c) Jacobi ensembles. The values of 
parameters considered are indicated in the caption. In fig. [31 we consider the density of transmission eigenvalues 
governed by JPD in eq. p3)) . Figures 4 (a), (b), (c) respectively show the density of transmission eigenvalues Tj, 
distribution of the dimensionless Landauer conductance G = J2j > average conductance (G) for several Af, N 
values. The analytical results for comparison in figs.[3l|4]are taken from [13, H^. Finally, in fig.[5]we show the density 
of eigenvalues for the /3-Wishart ensemble defined by eq. for M = 2N = 100, 6=1. The simulation results have 
been compared with the large N result of [ll| . 

The time for these simulations typically varied from 5 minutes for ^ 5 to 45 minutes for ^ 50 on a workstation 
with Intel Core2 6300 @1.86 GHz processor. The exact time, however, depends on the specific ensemble and the choices 
for the final convergence factor {F) [l^, the E and x windows, number of bins {ue, n^), and the production run steps 
(sp). For example, for fig. 3(a), (x = x), it took about 20 minutes for the choices F = 10~^, 1F£;=[0,500], ue = 500, 
Wx=[-7,7],nx = 70, sp = 5 X 10®. For fig. 5, (A = fy"^ = Px), the time taken was about 30 minutes for F = 10"^, 
We = [0,20000], ue = 400, Wy = [0,45], = 75, sp = 5 x 10^. 
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VI. CONCLUSION 



We proposed a method based on the Wang-Landau algorithm for numerical generation of the eigenvalue densities 
and related obscrvables of random matrix ensembles. This scheme produces the results for multiple /3 values in a 
single simulation, and is therefore advantageous compared to the conventional Monte-Carlo scheme. This feature 
also makes it ideal for investigating the /3-ensembles which have been of considerable interest in recent times. We 
demonstrated the utility of this method using several examples. 
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